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ABSTRACT 

We present a method for including steady-state gas flows in the plasma physics code Cloudy, which was 
previously restricted to modeling static configurations. The numerical algorithms are described in detail, to- 
gether with an example application to plane-parallel ionization-bounded H n regions. As well as providing the 
foundation for future applications to more complex flows, we find the following specific results regarding the 
effect of advection upon ionization fronts in H n regions: 

1. Significant direct effects of advection on the global emission properties occur only when the ionization 
parameter is lower than is typical for Hn regions. For higher ionization parameters, advective effects are 
indirect and largely confined to the immediate vicinity of the ionization front. 

2 The overheating of partially ionized gas in the front is not large, even for supersonic (R-type) fronts. For 
subsonic (D-type) fronts we do not find the temperature spike that has been previously claimed. 

3. The most significant morphological signature of advective fronts is an electron density spike that occurs at 
the ionization front whenever the relative velocity between the ionized gas and the front exceeds about one half 
the ionized isothermal sound speed. Observational evidence for such a spike is found in [Nn] /I6584 A images 
of the Orion bar. 

4. Plane-parallel, weak-D fronts are found to show at best a shallow correlation between mean velocity 
and ionization potential for optical emission lines even when the flow velocity closely approaches the ionized 
sound speed. Steep gradients in velocity versus ionization, such as those observed in the Orion nebula, seem to 
require transonic flows, which are only possible in a diverging geometry when radiation forces are included. 
Subject headings: gas dynamics, H II regions, numerical methods 



1. INTRODUCTION 

The classic early work on the effects of dynamics on the 
em ission structure of H n regions was carried out by H arring- 
ton ( |1977 ), who studied weak-D fronts in which the gas mo- 
tions are always subsonic with respect to the front. Within the 
fully ionized interior of the H n region the gas was found to be 
close to thermal and ionization equilibrium. Significant non- 
equilibrium effects induced by the dynamics are confined to 
the edge of the region, near the ionization front, where there 
exist large gradients in the radiation-field intensity and in the 
physical cond itions of the gas such as temperature and degree 
of ionization, ^iarrington found that the dynamics only had a 
small effect on the integrated forbidden line spectrum of the 
models he considered. However, there are various reasons to 
revisit such calculations now. 

First, ionization fronts are now studied in a diverse range 
of astrophysical contexts in which the clas sical, spheri cally 
symmetric, subsonic expansion studied by Harrington may 



ized regions, ranging in scale from cq metary knots in plan- 
etary nebula ( |L6pez-Martm et aT 2001 ) and phot oevaporated 



circumstellar disks in Hn regions (O'Dell 2001 and refer 



ences therei n) up to champagne flows in giant extragalactic 



H n regions dScowen et al.| |1998|) and the photoey aporation 



be the exception rather than the rule. Transonic photoevap- 
oration flows seem to be a ubiquitous feature of photoion- 
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of cosmological mini halos ( Shapiro & Raga 2001 ). In such 
flows, non-equilibrium effects will be somewhat more impor- 
tant than in s ubsonic weak-D fronts d ue to the higher veloci- 
ties involved. Sankrit & Hestei (2000) made a first attempt at 
detailed modeling of the emission structure of the flow from 
the head of the columns in Ml 6, using static equilibrium mod- 
els for the ionization structure. 

Second, continuous improvement in the spatial resolution 
and wavelength coverage of observations, together with ad- 
vances in theoretical and observational atomic physics, now 
allow a much more detailed comparison between model pre- 
dictions and spatially resolved observations of a multitude of 
emitting species. In this context, even moderate and localized 
changes in the predicted spectrum due to dynamical effects 
can be important. 

Third, a dynamical treatment allows the self-consistent cal- 
culation of the velocity field, which allows comparison with 
high-resolution spectral line profile observations that provide 
further constraints on the models. Furthermore, it permits 
a unified treatment of the entire flow from cold, molecular 
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gas, through the photon-dominated or photodissociation re- 
gion (PDR), and into the ionized re gion. Previou s studies 
of non-equ ilibrium models of PDRs ([London! Il978|; N atta & 
Hollenbach Tl998| ; |Stoerzer & Hollenbach||1998D have tended 
to treat the PDR i n isolation withou t consi dering the Hn re- 
gion in any detail. Richling & Yorke ( 2000 ) presented numer- 
ical radiation-hydrodynamic simulations of the photoevapo- 
ration of proplyd disks but the physics of the PDR was calcu- 
lated in a simplified manner. 

In common with most other photo ionization codes, Cloudy 
( Ferland et al. 1998 ; ^erlandj 200C ) has traditionally calcu- 
lated static equilibrium models in which time-dependent ef- 
fects are neglected, such as isochoric (constant density) or 
isobaric (constant pressure) configurations. The task of com- 
bining hydrodynamics with detailed simulation of the plasma 
microphysics can be approached from one of two angles. 
One method would be to add the atomic physics and radia- 
tive transfer processes to an existing time-dependent hydro- 
dynamic code. The other, which is the method pursued in this 
work, is to add steady-state dynamics to an existing plasma 
physics code. 

The current paper is the first in a series of three that will 
present detailed results from our program to include a self- 
consistent treatment of steady-state advection in a realistic 
plasma physics code. This first paper introduces the method- 
ology employed in the series as a whole and then concen- 
trates on the restricted problem of "weak" ionization fronts 
in a plane-parallel geometry. The second paper of the series 
includes the molecular reaction networks necessary for mod- 
eling the neutral/molecular PDR, while the third paper con- 
siders "strong D" ionization fronts, where the gas accelerates 
through a sonic point, as found in divergent geometries such 
as the photoevaporation of globules. 

We first discuss the general problem of advection in ion- 
ization fronts (section [|) . We then describe the modifications 
that have been made to the Cloudy photoionization code in 
order to treat steady-state flow (section |^). Results from a 
small sample of representative models are presented in sec- 
tion I] and the application of our results to observations of the 
Orion nebula is discussed in section g. Further technical de- 
tails of the physical processes and computational algorithms 
are presented in a series of appendices. 

2. ADVECTION IN IONIZATION FRONTS 

The classification of ionization fronts depends on the be- 
havior of the gas velocity in the frame of reference in whic h 



Williams & Dyson 1996 ). If the gas is magnetized, then the 
classification becomes more complicated since there are now 
three wave speeds to take into account (Alfven speed plus 
fast and slow magnetosonic speeds) instead of just the sound 



the ionization front is fixed (Kahn 1954; Goldsworthy 1961 ). 
In this frame, the gas flows from the neutral side of the front 
(denoted upstream) towards the ionized side (denoted down- 
stream). If the upstream gas velocity on the far neutral side is 
subsonic with respect to the front then the front is said to be 
D-type, while if it is supersonic the front is said to be R-type. 
A further distinction is made between those fronts that con- 
tain an internal sonic point, which are said to be strong, and 
those that do not, which are said to be weak. For example, 
a weak-R front will have supersonic velocities throughout the 
front, whereas in a strong-D front the gas starts at subsonic ve- 
locities on the neutral side, accelerating through the front to 
reach a supersonic exhaust velocity on the ionized side. When 
the downstream gas velocity is exactly sonic the front is said 
to be critical. There is also the possibility of a recombination 
front, in which the sense of the gas velocity is reversed and the 
flow is from the ionized side towar ds the neutral side, with a 
similar range of possible structures ( Newman & Axford|l968 



fast-mode weak-R front, etc ( 


Redman et al. 


1998; 


et al. 2000 


; Williams & Dyson 


2001). 



Williams 



These classification schemes were developed for plane- 
parallel fronts but will be approximately valid so long as the 
radius of curvature of the front greatly exceeds its thickness. 
This is usually the case since ionization fronts are in gen- 
eral very thin compared with the sizes of H n regions unless 
the ionization parameter is small (see below). What type of 
ionization front actually obtains in a given situation depends 
on the upstream and downstream boundary conditions of the 
front, in particular the upstream gas density and the down- 
stream gas pressure and ionizing radiation field, together with 
the large-scale geometry of the flow, which need not be plane- 
parallel. 

Since the gas velocity through the front is high for an R- 
type front, so is the flux of neutral particles that must be ion- 
ized for a given upstream density, which in turn requires a 
high ionizing flux at the downstream boundary. As a result, 
R-type fronts are usually transient phenomena accompanying 
temporal increases in the ionizing flux, such as the "turning 
on" of an ionizing source. In the most common case, the front 
will be propagating rapidly through slowly moving gas. In 
the limit of an extreme weak-R front, the gas velocity in the 
ionization front frame and density are constant throughout the 
flow. 

D-type fronts are more common and the limit of an ex- 
treme weak-D front corresponds to a static constant pressure 
front in ionization equilibrium. Weak-D fronts require a high 
downstream pressure and therefore are likely to be found in 
cases where the ionization front envelops the ionizing source. 
Strong-D and D-critical fronts, on the other hand, are con- 
sistent with the free escape of the downstream gas and hence 
apply to divergent photoevaporation flows, for example, from 
globules. 

Advection of material through the ionization front may be 
expected to have various effects on the emission properties of 
a nebula. In order to simplify the discussion, we will consider 
a plane-parallel nebula, illuminated at one face by a given ra- 
diation field and in a frame of reference in which the ioniza- 
tion front is at rest. This is illustrated in Figure [jjj. The gas is 
supposed to enter the front from the neutral side with velocity 
v n and to leave on the ionized side with velo city V[. Results 



from this simplified model are described in § |2.2| but we first 



discuss the relation between this model and real H n regions. 

2.1. Physical Context of Advective Fronts 

In this section, we consider two typical scenarios in which 
advective ionization fronts may be encountered and investi- 
gate to what extent they may approximated as steady flows in 
the frame of reference of the ionization front. In this discus- 



sion we follow |Shu| ( |1992| ) in denoting gas velocities in the 
frame of reference of the ionizing star by u, gas velocities in 
the frame of reference of the ionization front by u, and pattern 
speeds of ionization and shock fronts by U. 

2.1.1. Classical Stromgren Sphere 

The evolution of a classical Stromgren-type H n region in a 
const ant density medium has been descr i bed by ma ny authors 
(e.g., |Goldsworthyj|l958t |Spitzerj|l978|; ^hu||l992|; Dyson & 
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(c) 



Fig. 1. — Advective fronts: (o) Plane-parallel idealization; (b) Classical 
Stromgren sphere; (c) Photoevaporation flow. 



Williams 1997). If the ionizing source turns on instanta- 



neously, then the ionization front is initially R-type and prop- 
agates supersonically through the surrounding gas with little 
accompanying gas motion. By the time the ionization front 
reaches the initial Stromgren radius (where the rate of recom- 
binations in the ionized region approximately balances the 
ionizing luminosity of the source), the front propagation ve- 
locity has slowed to the order of the sound speed in the ionized 
gas and the front becomes D-type, preceded by a shock that 
accelerates and compresses the neutral gas. The initial R-type 
phase is very short (of order the recombination timescale), 
and consequently is of little observational importance. The 
structure of the region during its subsequent D-type evolution 
is shown in Figure [jj?. 

The propagation speed of the shock is roughly equal to that 
of the ionization front and also to the velocity of the gas in the 
neutral shell: 

U t - Uu - U D , (1) 

whereas the ionized gas expands homologously with a veloc- 
ity that increases linearly with radius, reaching half the speed 
of the ionization front at the front itself: 



t/if. 



(2) 



Hence, the velocity of the ionized gas immediately down- 
stream of the ionization front in the ionization front frame 
is 

v i = u i (R ] f)-U i[ = -0.5U il . (3) 

The ionization front propagates very slightly faster than the 
gas in the neutral shell, giving a small upstream neutral gas 
velocity in the ionization front frame of 

v n = u n - Uif < -2c 2 Jc{ ~ 0.1 kms -1 , (4) 

where c n , c\ are the isothermal sound speeds in the neutral and 
ionized gas, respectively. The evolution of the ionization front 
propagation speed can be described in terms of its radius as 

U * = ~ ' < 5 ) 



[4 0Rif/fli„it) 3/2 -l] 



where /?; n j t is the initial Stromgren radius. The Mach num- 
ber reached by the ionized gas just inside the ionization front, 
measured in the frame of reference in which the front is sta- 
tionary, is given by M — 0.5t/;f/ci. This is plotted in Figure || 
as a function of ionization front radius. During the lifetime 
of a typical O star (and assuming an ambient density of or- 
der 1 cm 3 ), the radius of a classical Hn region will expand 
by roughly a factor of 4, so, as can be seen from Figure [| 
M = 0.3-0.5 is typical of the majority of the evolutionary 
lifetime. 
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Fig. 2. — Evolution of the ionization front radius and the "exhaust" Mach 
number of newly ionized gas in the frame of reference of the advancing ion- 
ization front in a classical Stromgren sphere. These quantities are plotted as a 
function of time since the front reached the initial Stromgren radius, in units 
of the sound crossing time at that radius. 



2.1.2. Photoevaporation Flow 

Photoevaporation flows are very common in ionized re- 
gions, occurring whenever the ionization front is convex from 
the point of view of the ionizing source, thus allowing the 
ionized gas to freely stream away from the front. Examples 
include bright-rimmed globules and proplyds in H n regions 
and cometary knots in planetary nebulae. On a larger scale, 
blister-type H n regions can also be considered photoevapora- 
tion flows ( [Bertoldi & Drainej|l996| ). 

Figure [l|c show the structure of an idealized photoevapora- 
tion flow, r oughly corresponding to th e equilibrium cometary 
globules of |Bertoldi & McKee (199C). The neutral gas is as- 
sumed to be approximately at rest with respect to the ion- 
izing source and the ionization front to be D-critical with 
V[ — it[ — C{, The ionization front eats slowly into the neu- 
tral gas with U\ = —On = c n /(2e?) <K c;, and the Mach number 
reached by the ionized gas at the front will be 1 . Outside the 
front, the ionized gas accelerates as an approximately isother- 
mal w ind, as observed around the Orion proplyds (H enney & 
Arthur |1998[ )' 



2.2. Direct and Indirect Effects of Advection 

In order to provide some physical insight into advective ion- 
ization fronts and to guide the interpretation of the numerical 
simulations, we have developed a simple analytic model for a 
plane-parallel weak-D ionization front, in which the gas tem- 
perature is assumed to be a prescribed monotonic function of 
the hydrogen ionization fraction, x. With this assumption, and 
considering mass and momentum conservation, it is possible 
to find algebraic solutions for the density, gas velocity, and 
sound speed as functions of x (see Appendix |a|). These solu- 
tions form a one-parameter family characterized by the maxi- 
mum Mach number of the gas in the rest-frame of the ioniza- 
tion front, reached asymptotically as x — > I. If we now take 
into consideration the ionization balance and radiative trans- 
fer, one can find a solution for x(z), the ionization fraction 
as a function of physical depth, by solving a pair of ordinary 
differential equations. 

Apart from the maximum Mach number, Ai m , the solutions 
are found to depend on two dimensionless param eters, £ a d 
and t„. The first of these, £ a d (defined in equation |A18[ ), is 
roughly the ratio of recombination length to mean-free-path 
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Fig. 3. — Dynamic slab solutions from the simplified model developed 
in Appendix N. (a) T» = 30, low ionization parameter, fat IF. (b) r, = 
3000, high ionization parameter, thin IF. Electron density (solid line) and 
gas velocity (dashed line) are shown in each case for 6 models with M m = 
0.0, 0.2, 0.5, 0.7, 0.9, 0.99 (corresponding respectively to successive curves 
from right to left). Electron density is normalized by the fiducial density 
n m , velocity is normalized by the maximum sound speed c m) -and distance is 
normalized by the static Stromgren depth zo — see Appendix W for details. 



of ionizing photons in a D-critical front and is not expected 
to vary greatly between H n regions, having a typical v alue o f 
£ad —10. The second parameter, t* (defined in equation A17 ), 
is roughly the ratio of the thickness of the fully ionized slab to 
the thickness of the ionization front and is proportional to the 
ionization parameter at the ionized face: t* ~ 10 6 T, where 
Y = Fo/(noc). The global importance of advection for the 
system as a whole can be characterized by the parameter /l a( j, 
defined as the ratio of the flux of hydrogen atoms through the 
ionization front to the flux of i onizin g p hoton s at the illumi- 
nated face of the slab (equations A10 and \19 ). For small val- 
ues of /t a( j, its value can be approximated as /l a( j ^ a dA1 m /r*. 
This can be understood as follows: the local effects of ad- 
vection at the ionization front itself are always substantial (so 
long as M m is not too small), being of order £; a &M m . On the 
other hand, the partially-ionized zone occupies only a small 
volume compared with the fully ionized gas unless the ion- 
ization parameter is small, so the global effects of advection 
are reduced by a factor of r». 

Detailed results are calculated in AppendixIXIfor both a low 
ionization parameter model (t* = 30) and a nigh ionization 
parameter model (r» — 3000). The structure of these models 
for different values of M m is shown in Figure [| In the low- 



T model, the direct global effects of advection are expected 
to be significant and indeed the thickness of the ionized slab 
is reduced by 50% as M m approaches unity. In the high-T 
model (more representative of typical H n regions), the direct 
global effects of advection are expected to be negligible since 
^ad - 3 x 10~ 3 vVf m . However, we find that even in this case 
the thickness of the ionized slab varies by about 5% as Ai m 
is varied between and 1 . This is due to a pronounced peak 
that develops in the electron density distribution at the ion- 
ization front for M m > 0.5 (for a static model the electron 
density declines monotonically through the front). We also 
find that the ionization front becomes substantially sharper as 
M m is increased, which is due to a decrease in the ioniza- 
tion fraction for a given value of the optical depth to ionizing 
radiation. Both these indirect effects of advection have sig- 
nificant effects o n the emissivity profiles of optical emission 
lines (see Figure A19 of Appendix Q), especially those such 
as [Oi]/I6300 A and [S n] ^6716 + 6731 A that form close 
to the ionization front. One can also calculate the spectral 
profi les o f emi ssion lines from the models, as shown in Fig- 
ures A2C and |A2l| of the Appendix. Again, it is lines that 



form in the partially ionized zone that show the most inter- 
esting behavior. These may show RMS line widths roughly 
equal to the sound speed and significant velocity offsets, both 
due to the gas acceleration in the ionization front. The de- 
rived widths are roughly four times the thermal width of lines 
emitted by light metals. 

Although this analytic model has provided insight into 
some of the effects of advection, it is obviously deficient in 
many respects. Many physical processes have been ignored 
and in particular the use of a fixed temperature profile T(x) 
does not allow for the fact that T(x) itself may be affected by 
the advected flow. 

3. ADDING DYNAMICS TO CLOUDY 

In order to study the structure of advective ionization fronts 
in greater detail, we need to include a wide range of additional 
physics. This could be done in a variety of ways, for example 
by integrating through the steady state equations or by includ- 
ing source terms in a time-dependent hydrodynamic simula- 
tion. In each case, an implicit treatment of the source terms 
is necessary, as the many physical processes with timescales 
shorter than the dynamical timescale lead to the problem be- 
ing very stiff. 

In the present work, we have chosen to adapt the photoion- 
ization code, Cloudy, which already includes a comprehensive 
treatment of the physical source terms. Cloudy, in common 
with other traditional plasma codes, searches for equilibrium 
solutions of the ionization equations. In order to treat steady 
advective flows, we have included advective source and sink 
terms in the equilibrium balance equations. The effect of this 
is that the equilibrium search phase now in fact determines 
the implicit solution of the advective equations, and so treats 
short-timescale processes in a stable manner. 

In this section, we present the basic equations and outline 
the methods which we use to solve them. 

3.1. Equations 

Cloudy takes into account the conservation equations for 
each species and also the heating and cooling balance under 
the simplifying assumptions of constant density or constant 
pressure. This procedure can be summarized as balancing 
source and sink terms for the ionization and energy equations. 
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For the ionization equation in the static case this can be ex- 
pressed as 



dm 
dt 



0, 



(6) 



where dni/dt is the rate of change of the volume density of 
a particular ionization state, which in equilibrium is equal to 
zero. The R^j are the rates for ionization (where j is a higher 
state than /) and recombination (where j is a lower state than 
i). G, and 5, cater for processes not included within the ion- 
ization ladder, and are respectively the source of ions from 
such processes and the sink rate into them. A detailed discus- 
sion of the solution method for the ionization networks in the 
equilibrium case is given in Appendix 

The general Cloudy solution method works by a series of 
nested iterations. The innermost loop is the ionization net- 
work, external to this is the electron density iteration, which 
enforces charge neutrality, then the temperature loop, which 
enforces thermal balance. Finally, an optional outermost it- 
eration loop varies the density to achieve pressure (or more 
generally momentum flux) balance. The whole system is iter- 
ated until convergence within a given tolerance. 

Once dynamics is included, the continuity and momentum 
equations must be added to the set of equations to be solved, 
kinetic and internal energy transport and pressure work must 
be taken into account, and advection terms must be added to 
the ionization balance equations. For example, for a plane- 
parallel steady-state flow (the simplest case, but one that is 
applicable to blister Hn regions), the equations to solve in 
flux conservative form are 



7r<P«) = o, 

Ox 

—(p+pu ) -pa, 
ox 



d_ 

dx 



pu\w + —u 



= H-C, 



dx 



(n { u) =Gi + 2^ Rj^itij -tij S i + 2^ Ri^ j 



j*' 



(7) 
(8) 
(9) 

(10) 



Here, a is an acceleration e.g., gravity or radiation driving, w 
is the specific enthalpy w - s + p/p, where s is the specific 
internal energy, and H - C is heating minus cooling. Here 
the specific internal energy includes only the thermal energy 
of translation, so s = 3/2{p/p), as transfers from other physi- 
cal energy components (ionization energy, binding energy, vi- 
bration and rotation energy of molecules, etc.) are treated as 
explicit heating and cooling terms in the underlying thermal 
balance scheme. 

The advection terms have the general form of V ■ (n,y) (for 
steady state), where v is the advection velocity. This can be 
written as 

V ■(n i v) = nv-V(ni/n). (11) 

3.2. Differencing Scheme 

Although it is possible to solve the ionization equations in 
an explicitly time-dependent way, this is not the best way to 
proceed. The photoionization terms in the steady-state solu- 
tion will often have very short timescales, so stability con- 
straints would limit the time step to this short photoionization 
timescale, and hence cause an extremely slow convergence of 



the iterative scheme. Instead, we take advantage of the cur- 
rent algorithm used by Cloudy and difference the equations 
implicitly. Such an implicit scheme has the advantage that the 
timestep is not limited to the shortest ionization or recombina- 
tion time, which is clearly unsatisfactory for an astrophysical 
system in which the dynamical timescales are general much 
longer than the ionization or recombination timescales. 

At iteration m, the advection terms may be approximated 
based on the value in the present zone and an upstream value 
in the previous full iteration as 

d x^(z)-xf- l (z-Az) 

dz Az 
where x™(z) is the value of x-, = n;jn at position z at the mth 
iteration of the scheme, and Az is an adjustable advection 
length. For the first iteration no upstream values are avail- 
able so no advection terms are included in the equations. It is 
useful to define the look-back operator 

L Az [x?(z)]=x?-\z-Az) (13) 

for values, such as x%, given per unit material. Values speci- 
fied per unit volume need to be scaled to a conserved variable 
before the look back is applied, so that 

L Az [n?(x)]=n m x?-\z-Az). (14) 

This may be thought of as a first-order Lagrange-remap solu- 
tion for the advection equation. 

For the scheme discussed in Appendix for the ionization 
ladder, the advective terms may then be included simply as an 
additional source term, 

d = nv^-^ , (15) 



and sink rate 



Az 



Az 



(16) 



in the linearized form of the equations, and iterated to find the 
non-linear solution as in the time-steady case. 

The energy balance equation may be treated in a similar 
manner. Using the mass conservation equation (eq. [^), we 
have from equation (p|) 



r, 1 ' 

pv.W I u; + —v* 

which may be differenced as 

(w + jv 2 ) - L A: ( 



H-C, 



w + 



pv- 



Az 



= H-C. 



(17) 



(18) 



The terms on the left side of this equation may then be treated 
as additional heating and cooling terms in the temperature 
solver. 

The continuity and momentum equations are more easily 
dealt with. The continuity equation is taken into considera- 
tion simply by eliminating the velocity v in terms of p in all 
equations using the substitution 

pvr d = constant, (19) 

which comes from integrating the general form of equation 
where d — 0,1,2 indicates plane parallel, cylindrical or spher- 
ical geometry, respectively. The initial condition is given at 
the illuminated face. 

The dynamical pressure, which appears in the momentum 
equation, is taken into account by adding the ram pressure 
term pv 2 to the total pressure. 



0.1 



0.01 



0.001 



10-" 



10" 



"i 1 r 



~ 1 Akd/lO^cm - 
ei - 

f -2 ■ ■ 



10 



20 



30 40 50 60 
Iteration Number 



70 



SI) 



DO 



Fig. 4. — Convergence behavior of an example model as a function of iter- 
ation number, m. Solid line: advection length Az (units of 10 15 cm); dashed 
line: convergence error, ei; dotted line: discretization error, £2. 



Varying the Advection Length — The advection length, Az, in 
this scheme determines the manner in which different pro- 
cesses are treated. The differencing we have chosen has the 
effect that processes far more rapid than Az/v are treated as 
in static equilibrium, while slower processes are followed ex- 
actly. The correct steady-state solution is found in the limit 
Az — > 0, but the smaller the advection length chosen, the 
longer the system will take to reach an equilibrium state. The 
natural procedure is then to use a first iteration solution by 
ignoring the advection terms, and then gradually decrease Az 
until Az ~ Azgrid, where Az gr id is the size of a spatial zone in 
the simulation, at which stage the treatment of advection will 
be as accurate as that of photoionization. 

In order to track the convergence of the models and to deter- 
mine when to reduce the advective timestep, we monitor the 
behavior of two error norms, e\ and ex. Both are calculated 
as the squared norm over all zones, z, and ionic/molecular 
species, i. The first of these norms is the convergence error, 
defined as 



Az/v 



(20) 



which measures the difference between the model solutions 
for the last two iterations. The second is the discretization 
error, defined as 



£2 = 



Bj - L Az (nj) Hi - L(Az/2)("i) 



Az/v 



Az/2v 



(21) 



which measures the accuracy of the present estimate of the 
advective gradients in the solution compared to an estimate 
with half the advection length. 

If £\ <k 62, then the solution is converged with the present 
Az, while if the errors in the Lagrangian estimate of the gradi- 
ent of the value are still significant, then the timestep should 
be decreased. Cutting Az when < O.le? produces substan- 
tial improvements in the rate of convergence of the advective 
solutions. However, it can still take a substantial number of 
iterations to reach equilibrium for large advection velocities. 

An example of the way in which Az, & , and ei vary during 
a model calculation is given in Figure |] (which corresponds 
to model ZL009 discussed in the following section). It can 
be seen that while Az at j remains constant, the convergence er- 
ror, ei , decreases with each iteration, while €2 hardly changes 
after the first iteration for a given Az. When e\ / £2 falls to a 



low enough value, then all physical processes that occur on 
timescales longer than At = Az/v have converged, so the ad- 
vection length can be reduced. This has the effect of lower- 
ing €2 but also temporarily increases e\ due to the release of 
shorter-timescale processes from strict local equilibrium, so 
several iterations must be carried out at the new value of Az. 
This procedure is continued until €2 has fallen to a sufficiently 
low value. 

4. RESULTS 

In this section we present results for selected advective 
ionization fronts calculated using Cloudy, all using a plane- 
parallel slab geometry. The principal input parameters for the 
models are the hydrogen number density, no, gas velocity, uq, 
hydrogen-ionizing photon flux, Fq, all specified at the illumi- 
nated face. The spectral distribution of the incident radiation 
field was assumed to be a black body with effective tempera- 
ture, r e ff- All these parameters are shown in Table [j] for the 
three models presented here. The gas-phase elemental abun- 
dances for all the models were set at the standard ISM val- 



ues ( [Baldwin et al.||1996|) and Orion-type si licate and graphite 
grains were included ( paldwin et al. 1991). Since this paper 
is concerned with the effects of advection on the ionization 
front, all molecules were turned off and the integration was 
stopped when the electron fraction fell below 10~ 3 . The in- 
clusion of molecular processes in the PDR will be described 
in a following paper. The models also include an approxi- 
mate treatment of a tangled magnetic field (see Appendix^), 
characterized by the field strength at the illuminated face, Z?q. 

One further physical process was disabled in these mod- 
els: the radiative force due to the absorption of stellar contin- 
uum radiation (principally by dust grains). This was done for 
purely pragmatic reasons, since the inclusion of this process 
for high ionization parameter models makes it very difficult to 
set the approximate desired conditions at the ionization front 
by varying conditions at the illuminated face. Models that do 
include this process are discussed further in section p3[ 

The first two models, ZL009 and ZH007, are weak-D fronts 
with low and high ionization parameter, respectively, with pa- 
rameters similar to the toy models discussed in Appendix^. 
The velocity at the ionized face, uq, was set somewhat below 
the isothermal sound speed in order to avoid the possibility 
of gas passing through a sonic point during an intermediate 
iteration (transonic fronts will be considered in a following 
paper). The third model, ZH050, is a weak-R front in which 
the gas velocity relative to the front is supersonic throughout. 
Such R-type fronts are likely to be transient and thus of lim- 
ited observational significance. Nevertheless, this model is 
useful since it provides a stringent test of our simulations in 
the limit of high advective velocities. 

As well as the advective models, we also calculate equiv- 
alent static models for each of the three cases considered, 
which have constant pressure for comparison with the weak- 
D fronts or constant density for comparison with the weak-R 
front. 

4.1. Low Ionization Parameter, Weak-D 

This model, ZL009, has physical parameters that are in- 
spired by those of cometary globules in planetary nebulae 
such as the Helix, although the geometry is plane-parallel 
rather than spherically divergent. The ionization parameter of 
the model is very low (T = 3.33 x 1CT 5 ), which accentuates 
the global effects of advection and also leads to the ionization 
front thickness being comparable to the Stromgren thickness 
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TABLE 1 
Model Parameters 



Parameter 


ZL009 


ZH007 


ZH050 


log(n /crrr 3 ) 


3.5 


4.0 


4.0 


Ho/krn s 


-9.0 


-7.0 


-50.0 


log(F /cm" 2 ) 


9.5 


13.0 


13.0 


log(r cff /K) 


4.6 


5.0 


4.6 


log(.Bo/Gauss) 


none 


-4.0 


-4.0 


log T 


-4.5 


-1.5 


-1.5 


r» 


9.9 


6800 


6800 


Mm 


0.84 


0.73 


4.01 


6id 


8.1 


13.8 


13.8 


<*ad 


2.2 


0.0015 


0.008 




■ T <: (advective 

■ T c (static) 




l-lO" 2.10 1 " 310 is 4-10" 

Depth, 2 (cm) 

Fig. 5. — Structure of model ZL009 as a function of depth from the illu- 
minated face, (a) Velocity and isothermal sound speed, (b) Gas temperature, 
(c) Number density of hydrogen nucleons. (d) Hydrogen ionization fraction, 
(t?) Helium ionization fractions. Panels b—e show the advective model results 
(solid line) and the results from an equivalent static constant pressure model 
(dashed line). 
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Fig. 6. — Emissivity structure of model ZL009. As Figure ^| but showing 
volume emissivity of (a) Ha ,16563 A, (b) [Nil] -16584 A, (c) [Oi] -16300 A, 
(d) [S II] -16716 + 6731 A, and (e) [Om] -15007 A. 
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Fig. 7. — Face-on emission line profiles of model ZL009 including ther- 
mal broadening: Ha -16563 A (thick solid line), [Nil] -16584 A (thin solid 
line), [Oi] -16300 A thick dashed line), [S II] -16731 A (thin dashed line), and 
[O III] -15007 A (thick dotted line). All lines are normalized to their peak in- 
tensities. 



of the ionized layer. 

As can be seen from Figure |[ the advection has large ef- 
fects on the model structure as would be expected given the 
large value of /t ac j. The depth of the ionized region is reduced 
by more than a factor of two with respect to the static model, 
with concomitant reductions in the brightness of the hydro- 
gen recombination lines (see Table The collisional lines 
of [Oi], [Sii], and [Nn] are also reduced in intensity, albeit 
to a lesser degree ([Om] emission from this model is neg- 
ligible due to the low ionization parameter). Interestingly, 



the intrinsic Balmer decrement is increased in the advective 
model, which gives Ha/H/3 = 3.29 (reddening by internal 
dust is not included in the line ratios given in Table This is 
because the temperature at low ionization fractions is signif- 
icantly higher in the advective model, which leads to a non- 
negligible collisional contribution to the Balmer emission that 
preferentially excites Ha. The temperature in the more highly 
ionized zone is also increased somewhat by the effects of ad- 
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vection. 

Also listed in Table ^| are the mean velocities and full- 
width-half-maxima of the different emission lines, calculated 
assuming that the model is observed face on. The emission 
line profiles are illustrated in Figure ^. The contribution from 
each computational zone to the line profile is thermally broad- 
ened using the local temperature and the atomic weight of 
the emitting species. As a result, the Ho- line is significantly 
broader than the collisionaly excited metal lines. The [No] 
line is blue shifted with respect to [Sn] and [Oi], which is 
due to the acceleration of gas through the ionization front, as 
can be seen in Figure ^p. The Ha line has two components: a 
broad blue-shifted component due to emission from the ion- 
ized flow, and a narrow component at zero velocity, caused 
by a subsidiary peak in the electron density that occurs at low 
ionization fractions where the temperature changes sharply. 

For the [S n] line, the thermal broadening and the gas accel- 
eration both contribute in equal amounts to the predicted line 
width. For the lighter metal lines, the acceleration broadening 
remains roughly the same but the thermal width is increased 
somewhat. For Ha, the thermal broadening dominates. 

In this model, the gas velocity and Mach number increase 
monotonically as the gas flows from the neutral to the ionized 
side, reaching a maximum Mach number of Ai m - 0.84 at the 
illuminated face. As a result, emission lines from more highly 
ionized species are more blue-shifted (see Table ||). However, 
the effect is slight with only 2.5 km s velocity difference 
between [O i] and [O m] . 

4.2. High Ionization Parameter, Weak-D 

This model, ZH007, has physical parameters (see Table [[]) 
inspi red by the central region of the Orion nebula (B aldwin 
et al. 1991 ). The density is only slightly higher than in ZL009 
but the ionizing flux is much higher, resulting in a far higher 
ionization parameter (T = 3.3 x 1CT 2 ). The ionization front is 
much thinner than the depth of the fully ionized slab. The lo- 
cal advection parameter, ^ a( j (see Section || and Appendix |a|), 
is somewhat higher than in the previous model due to the 
hotter gas temperature and softer ionizing spectrum but the 
much higher value of t* that accompanies the higher ioniza- 
tion parameter means that the global advection parameter, /i a d, 
is very small. 

The structure of the advective model as a function of depth 
into the slab from the illuminated face is shown by solid lines 
in Figure ^| For some panels, an equivalent constant-pressure, 
static model is also shown (dashed lines). The indirect effects 
of advection are much greater than the direct loss of 0.1% 
of the incident ionizing flux due to the ionization of fresh gas. 
The largest effect is a roughly 3% increase in the mean density 
in the ionized zone (the densities at the illuminated faces are 
set equal in the static and advective models) due to the varying 
importance of ram pressure as the velocity varies through the 
slab. Due to the difference in density dependence of recom- 
bination and dust absorption, this leads to a slight decrease in 
the ionized column density together with an increase in the 
emission measure, caused by a reduction in the fraction of 
ionizing photons that are absorbed by dust grains. 

Figure || shows some of the same quantities as in Figure || 
but this time plotted against the electron fraction, jt e = n e /n. 
This effectively 'zooms in' on the ionization front transition 
itself, allowing one to appreciate details of the structure that 
are not apparent in the plots against depth. Note that the more 
neutral gas is on the left in Figure |[ whereas it was on the 
right in Figure [s| In Figure ^p, which shows the temperature 
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Fig. 8. — Structure of model ZHC1Q7 as a function of depth from the illumi- 
nated face. All panels as in Figure nl 



profiles, it can be seen that the initial heating of the gas as it 
is ionized is more gradual in the advective model than in the 
static model due to the photo-heating timescale being longer 
than the dynamic timescale. However, once the ionization 
fraction exceeds about 20% the photo-heating rate exceeds 
the static equilibrium model because the neutral fraction at 
a given value of the ionizing flux is higher than in the static 
case. This produces a characteristic overheating, which is of- 
ten seen in advective fronts. In the current case it is relatively 
modest, producing an extended 10 4 K plateau for x e = 0.3- 
0.9. For x e > 0.9 the gas approaches static thermal equi- 
librium again and the two curves coincide with temperature 
variations determined by radiation hardening and the varying 
importance of the diffuse field. Figure ^ shows the variation 
in the gas density and it can be seen that the density on the 
far neutral side is significantly higher in the advective model. 
This is a result of the lack of ram pressure on the neutral side 
(see Figure ^), which means the thermal pressure must in- 
crease to compensate. The magnetic field at the illuminated 
face in this model was set to be 0. 1 mG, implying a negligible 
contribution to the total pressure in the ionized gas. However, 
on the far neutral side, where the density is higher the field is 
much larger (assumed to grow with compression as B oc p 2 ^ 3 , 
see Appendix ^), so that the magnetic pressure becomes ap- 
preciable, limiting the density in the cold gas. The inferred B 
in the neutral gas is similar to ob served values in the neutral 



veil of Orion (Iroland et al. 1989) 



The line emissivities as a function of depth for the ZH007 
model are shown in Figure HQ. The main change between the 
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TABLE 2 

Emission Line Properties from Weak-D Photoionization Models 







Model ZL009 


Model ZH007 






Line Intensity 


Line Intensity 


Line 


A 


Static Advect v Av 


Static 


Advect v Av 


Ho- 


6563 


2.95 3.29 -5.89 21.6 


2.92 


2.92 -6.46 20.0 


[Ol] 


6300 


2.40 4.91 -5.04 6.3 


3.5 (-3) 2.8 (-3) -6.99 6.4 


[Nil] 


6584 


3.97 6.23 -6.07 6.9 


0.42 


0.43 -7.56 6.0 


[S II] 


6731 


3.52 5.33 -5.17 5.3 


0.04 


0.03 -7.29 4.3 


[Om] 


5007 


0.05 0.03 -7.50 6.1 


5.21 


5.21 -6.32 5.0 


US 


4861 


-2.868 -3.382 


0.435 


0.440 


Note. — 


Line intensities are all calculated for 


a face-on 


orientation and are 


given 


relative to H/3 except for HyS itself, which 


is log(intensity) in units of 


erg cm - s 










Fig. 9. — Structure of model ZH007 as a function of electron fraction, n e /n. 
Panels (a)-(c) as in Figure H. (d) Ionization fractions of H + (medium weight 
line), He + (thick line), and He 2+ (thin line), (e) Partial contributions to the 
total pressure: thermal gas pressure (medium weight line), magnetic pressure 
(thin line), and ram pressure (thick line). 



static and advective models is the sharp peak at the ionization 
front that is seen in the emissivity of lines from singly-ionized 
species such as [N ill. This peak is due to the electron density 
peak seen in Figure [| and discussed in Appendix |a[ 

The integrated emergent emission line spectrum (Table ||) 
is barely affected by the advection. As a result of the increase 
in emission measure discussed above, although the ionization 
front moves appreciably inward (Figure [|f), the Balmer line 
flux is actually slightly higher in the advective model than in 
the static model. On the other hand, lines that form in the 
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Fig. 10. — Emissivity structure of model ZH007 in units of 2 X 
10~ 16 erg cirT 3 s -1 ). Panels as in Figure ra 



ionization front itself, such as [S n] and [Oi], are somewhat 
reduced in intensity in the advective model due to the nar- 
rowing of the ionization front by advection (see discussion in 
Appendix 

Unlike in the the low ionization parameter model of the 
previous section, in this model the velocity and Mach num- 
ber do not increase monotonically from the neutral to the ion- 
ized side. Instead, the maximum Mach number (A4 m = 0.73), 
which also corresponds to the maximum velocity, occurs at 
the He°/He + front. Although this does not correspond to the 
maximum temperature, it is a maximum in the sound speed 
due to the decrease in the mean mass per particle when He 
is ionized. There is a second, slightly lower, maximum in 
the sound speed, Mach number, and velocity just inside the H 
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Fig. 1 1 . — Face-on emission line profiles of model ZH007 including ther- 
mal broadening. Hff/16563 A (solid lines), [Nll]/16584 A (dashed lines), 
[Ol]i6300 A (dotted lines), and [S II] ,16731 A (dot-dashed lines). Vertical 
axis units are arbitrary. 



ionization front where radiation hardening causes a tempera- 
ture maximum. The velocity also starts to increase again very 
close to the illuminated face. 

Due to this complex velocity structure, the trends of 
blueshift with ionization parameter are less clear in this 
model, as can be seen from Table || and from Figure [TT|, which 
shows simulated emission line profiles. 

4.3. High Ionization Parameter, Weak-R 

This model, ZH050, is identical to ZH007 except that the 
gas velocity at the illuminated face is set to 50 km s _1 , produc- 
ing a weak-R front. The model structure as a function of depth 
is shown in Figure 12 and as a function of electron fraction in 
Figure In both cases, the advective model is now com- 
pared with a constant-density static model as opposed to the 
constant-pressure model that was used for comparison with 
the weak-D models. As can be seen from Figure ^j\p and c 
the velocity and density are roughly constant across the front, 
as is expected in the weak-R case. The extremal Mach num- 
ber in the front (which for R-type fonts is a minimum, see 
Appendix|A|) is M m = 4.01, which occurs on the ionized side 
of the front (see Fig. file) at x e 0.93. 



4.4. Temperature Structure of the Ionization Fronts 

Figure [j~2|Z? shows that our weak-R model has a pronounced 
temperature spike at the ionization front, together with a 
smaller spike at the He ionization front. Figure pj| b shows 
that this is a more extreme manifestation of the overheating 
in the ionization front that was seen in the weak-D model 
(Fig. However, even in this case of a rapidly propagating 
front we do not find the overheating to be very great, reach- 
ing a maximum temperature of only 13, 000 K. Other authors 
have found more pronounced over-heating effects in dynamic 
ionization fronts (fcodriguez-Gaspar & Tenorio-Tagle 1998 



|Marten & Szczerba 1997|) but for diffe rent values of the phys 
ical parameters. Marten & Szczerbl] found temperatures as 
high as 20, 000 K behind R-type fronts moving at a substan- 
tial fraction of the speed of light, which would be difficult 
to model usi ng our steady-state code. R odriguez-Gaspar & 
Tenorio-Tagle studied the time-dependent evolution of an H n 
region after the turning-on of an O star, finding temperatures 
of order 15,000 K behind the ionization front soon after its 
transition from R to D-type. Their higher temperatures may 
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Fig. 12. — Structure of model ZHOSp as a function of depth from the illu- 
minated face. All panels as in Figure [ 



be due to these authors including fewer cooling processes than 
are inc luded in Cloudy (for further discussion, see W illiams & 
Dyson ^OOl , and referen ces therein). The greater width of the 
temperature spike in the Rodriguez-Gaspar & Tenorio-Tagle 
models is due to the fact that they were considering lower den- 
sities, so the gas moves further from the ionization front in the 
time it takes to cool to the equilibrium temperature. 



5. DISCUSSION 

In this section, we look for evidence of advective effects 
in one of the closest and best-studied H n regions, the Orion 
nebula. We concentrate on the clearest signatures of advection 
to emerge from our simulations: the electron density spike 
and the ionization-resolved kinematics. 

5.1. Electron Density Structure of the Orion Bar 

One firm prediction of the advective model that differs from 
the static case is the existence of a sharp electron density peak 
at the position of the ionization front. This peak manifests 
itself most clearly in the emission of the [Nn] /16584 A line 
(see Figure [l(]), producing a narrow emissivity peak at the 
edge of the broader peak that comes from the neutral Helium 
zone. In Figure Jl4] we show an [N n] image of the bright bar 
region in the Orion nebula that shows evidence of just such a 
structure. 2 The bar is believed to be a section of the principal 
ionization front in Orion that is seen almost edge-on. It can 

2 This image is based on data obtained with the WFPC2 instrument on the 
Hubble Space Telescope, provided by C. R. O'Dell. 
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Fig. 13. — Structure of model ZH050 as a function of electron fraction, 
n c /n. All panels as in Figure except for (c), which shows the isothermal 
Mach number. 



be seen as a diffuse strip of [Nn] emission (width 10" =s 

6 x 10 16 cm) stretching from the top-left to bottom right of 
the image. The principal source of ionizing radiation is the 

07 V star 6 l OriC, located off the image to the NW. A sharp, 
bright edge to the emission can be seen on its SE side along a 
considerable fraction of the bar, which may correspond to the 
electron density peak. 

The geometry of the nebula in the vicinity of the bar is 
far more complicated than the plane-parallel geometry as- 
sumed in the models, making a quantitative comparison dif- 
ficult. The bar probably consists of at least two overlapping 
folds in the ionization front and its appearance is also affected 
by protruding fingers of neutral gas and interactions with the 
HH 203/204 and HH528 jets. However, the straight region 
of the bar to the NE of the HH 203/204 bowshocks show a 
relatively simple structure, which we will attempt to compare 

i emission 



with our model predictions. We present in Figure 15 



line spatial intensity profiles along a short section of narrow 
slit parallel to the bar, with position as indicated in Figure [14l 
Comparison of Figure |5] with the model profiles of Figure]^ 
shows good general agreement. 

The electron density in the bar region has been mea- 
sured from t he [S nl /16716 + 6731 A line ratio to be around 



4000 cm- 3 ( |Wen & Q'Dell| |1995j parcia Diaz & Henney 
2003 ), whereas model ZH007 predicts a value of ^ 8000 crrT J 



for the [S n]-derived density. On the other hand, the incident 
ionizing flux can be estimated to be about 5 x 1 1 2 cirT 2 , which 
is half that of the model ZH007. Thus the model has the same 
ionization parameter as the bar and the results can be directly 



compared by multiplying the model lengths by a factor of two 
such that the spatially broad component to the [N n] emission 
is predicted to have a thickness of ^ 4 x 10 16 cm ^ 7", in rea- 
sonable agreement with what is observed. The electron den- 
sity peak at the ionization front is predicted to have a thick- 
ness of =! 2 x 10 15 cm =i 0.3", which is also very close to the 
observed thickness of the narrow ridge in [N n] (note that this 
thickness is fully resolved by the HST, which has an angular 
resolution of 0.1" at optical wavelengths). 

5.2. Velocity-Ionization Correlation in the Orion Nebula 

The mean velocity of different optical emission lines from 
the core of the Orion nebula has long been known t o cor- 
relate with the ionization potential of the parent ion ( Kaler 
19671; D'Dell et alj h993|foenney & Q'D e l]| |1999[ O 'Dell 



et al. [2004 P oi et alj[2004| ). Lines from" more highly ionized 
species such as [O in] /15007 A are blue-shifted by approxi- 
mately 10 km s with respect to the molecular gas of OMC- 
1, which lies behind the nebula, with intermediate-ionization 
species such as [N n] /16584 A being found at intermediate ve- 
locities. 

The results of section Q show that just such a correlation 
can be qualitatively reproduced by our models. However, on 
closer inspection, many significant differences are revealed 
between the model results and the Orion observations. A 
much clearer velocity-ionization correlation exists in model 
ZL009 (low ionization parameter) then in ZH007 (high ion- 
ization parameter, more pertinent to the Orion nebula), as can 
be seen from Table ^| and Figures |and|ll} This is not sur- 
prising since the emission in the high ionization parameter 
model is dominated by ionized equilibrium gas, where veloc- 
ity changes are rather small and are driven by variations in the 
thermal balance (the complex interplay of radiation hardening 
and the excitation of different coolant lines) rather than be- 
ing directly caused by ionization changes. Furthermore, both 
the magnitude of the velocity shifts seen in the models and 
the broadening induced by the gas acceleration are somewhat 
smaller than is observed in Orion. 

In order to better reproduce the observations, what is re- 
quired is a means of continuing the acceleration of the gas 
inside the body of the nebula, where hydrogen is fully ion- 
ized. There are two means by which this might be achieved: 
first, by including the continuum radiation force upon the ion- 
ized gas/dust mixture, or, second, by considering a transonic 
strong-D ionization front in a non-plane, divergent geometry. 

The results of a model that includes the continuum radia- 
tion force are shown in Figure |l6[ The thick line in the lower 
panel shows the integral along the radiation propagation di- 
rection of the radiative force per unit volume, f m &, that acts on 
the material in the flow. 3 In response to this radiative forcing, 
the total pressure increases with depth and, since the flow is 
subsonic, the gas pressure and density increase in the same 
direction. By mass conservation, this leads to an acceleration 
of t he g as towards the illuminated face, as can be seen in Fig- 
ure |i"6p (compare Figure ^Ja, where this process is absent). 

However, this acceleration of the fully-ionized gas does not 
have a very large effect on the mean velocity of the emis- 
sion lines, as can be seen from Figure |l7| This is because 
the higher-velocity gas represents only a small fraction of the 
total emission, even for the [Om] line. Although the model 

3 The dust-gas drift velocity is always much slower than the flow velocity 
in these models, so it is valid to suppose that the gas and dust dynamics are 
perfectly coupled. 




Fig. 14.— Negative HST WFPC2 image of the Orion bar in the [Nil] i6584 A line. Image dimensions are 150 X 120". 
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the graph at an offset of = -120". 
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Fig. 16. — Structure of model ZHR012, which includes the continuum radi- 
ation force, (a) Velocity and isothermal sound speed, {b) Partial contributions 
to the total pressure: thermal gas pressure (medium weight line), ram pres- 
sure (thin line), integrated radiative force (thick line), resonance line radiation 
pressure (dashed line), and magnetic pressure (dotted line). 



ZHR012 does show a clear relation between velocity and ion- than 1 km s 1 between [O i] and [O m], compared with an ob- 
ization (unlike ZH007), the gradient is very small, being less served difference of 5 to 10 km s _1 in the Orion nebula. Also, 
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Fig. 17. — Mean velocities of common optical emission lines from plane 
parallel advective models in face-on orientation. The ionization potential of 
the parent ion increases from left to right. 



any deviations from a strictly face-on orientation of the ob- 
server would tend to reduce the observed gradients. 

In order to reproduce the observed kinematics of the Orion 
nebula, one needs a strong acceleration of the gas in the region 
between the hydrogen and helium ionization fronts, since this 
is the region where the ionization of heavy elements such as 
oxygen and sulfur is changing most swiftly. We have shown 
that it is not possible to achieve this with a plane-parallel 
model. In such a geometry, gas acceleration requires either 
a gradient in the sound speed or the application of a body 
force, neither of which are present with the required mag- 
nitude in the relevant region. Strong changes in the sound 
speed only occur at greater depths, in the hydrogen ionization 
front, while the effective gravity associated with continuum 
radiation pressure acts mainly at shallower depths, where the 
heavy-element ionization is not changing. 

On the other hand, as we will show in a following paper, 
the requisite acceleration is a natural consequence of mod- 
els in which the flow is divergent (either spherical or cylin- 
drical) rather than plane-parallel. Although such models are 
obvio usly relevant to such objects as proplvds ( Henney & 
Q' DellW99|) and photoevaporating globules (B ertoldi & Mc- 
Kee |1990[ ), it is not so apparent that they should apply to the 
large-scale emission from the Orion nebula, which has been 
traditionally visualized as a bowl-like cavity on the near side 
of the molecular cloud OMC-1. However, three-dimensional 
recons truction of the shape of the ionization front ( Wen & 



0'Delf[n?95 ) indicates that the radius of curvature of the front 
is larger than its distance from the ionizing star. In such a case, 
the divergence of the r adiation field c an lead to a (weaker) di- 
vergence of the flow (Henney 2003). Another possibility is 
that the mean flow is the superposition of multiple divergent 
flows from the many bar-like featur es that have been found in 
the nebula flO'Dell & Yusef-Zadeh||2000| ). 



6. CONCLUSIONS 

In this paper, we have developed a method for including 
the effects of steady mater ial flows in the plasma physics 
code Cloudy ( Ferland 200C ), which was previously capable 
of modeling only static configurations. We have presented a 
detailed description of the numerical algorithms and an ex- 
ample application to the restricted problem of plane-parallel 
ionization-bounded H n regions where the flow does not pass 
through a sonic point (weak fronts). The most important con- 
clusions from our study are as follows: 



1. The local effects of advection are most important in 
those regions of the flow where physical conditions are 
strongly varying over short distances, such as in the hy- 
drogen ionization front. 

2. The global effects of advection on an Hn region de- 
pend on the relative thicknesses of the ionization front 
and the region as a whole, which is a function of the 
ionization parameter. Only for low values of the ion- 
ization parameter, such as are found in cometary knots 
of planetary nebulae, do we find a significant direct ef- 
fect of advection on the emission properties integrated 
over the entire the region. For higher ionization param- 
eters, more typical of H n regions around O stars, the 
effects of advection are indirect and more localized. 

3. One such indirect effect is a modification of the temper- 
ature structure in the ionization front due to the over- 
heating of partially ionized gas. However, we find 
the magnitude of this ef fect to be much le ss than has 
previously been claimed (|Osterbrock||l989|; R odriguez- 
Gaspar & Tenorio-Tagle |199*8 ), probably due to our 
more realistic treatment of heating and cooling pro- 
cesses. For weak D-type fronts (subsonic flows), the 
temperature reached in the front does not exceed the 
equilibrium temperature in the fully ionized gas (see 
Figure ^>). Even for supersonic R-type fronts, the peak 
temperature is only about 20% higher than the equilib- 
rium ionized value (see Figure |l3|fo). The temperature 
increase causes an increase in the peak emissivity of 
the [O i] /16300 A line, but the total emission of this line 
tends to be less than in an equivalent static model be- 
cause advection acts to sharpen the ionization front and 
hence decreases the width of the zone where [Oi] is 
emitted. 



Another indirect effect of advection in high ionization 
parameter regions is the production of a sharp spike in 
the electron density, which occurs at the ionization front 
whenever the peak Mach number in the flow exceeds 
about 0.5. This spike does not occur in static models 
and its existence can be shown analytically to be a con- 
sequence of the exchange between thermal pressure and 
ram pressure as the g as is a ccelerated through the front 
(Appendix |a|, Figure A18). As such, it can serve as a 
useful diagnostic for the presence of advection in ion- 
ization fronts, best observed in the [N n] /16584 A line, 
for which advective models predict a two-component 
structure (Figure |To|?7): a broad peak of emission from 
the equilibrium H + -He zone plus a narrower peak from 
the « e spike at the ionization front. Just such a structure 
is seen in HST images of the Orion bar (Section 5.1), 



suggesting that advection is important in that region. 

Finally, the advective models provide a mapping be- 
tween the ionization state of the gas and its velocity 
and can hence be used to predict kinematic profiles of 
different emission lines, which can be compared with 
spectroscopic observations. We make such a compar- 
ison in the case of the Orion nebula but find that the 
plane-parallel models presented in this paper are utterly 
inca pabl e of explaining the observed kinematics (Sec- 
tion 5.2). The observations seem to demand a strong 



acceleration of the gas in the region between the hydro- 
gen and helium ionization fronts, whereas the only ac- 
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celeration mechanisms that can act in the plane-parallel 
models either occur at greater depths (gas heating in 
the hydrogen ionization front) or at shallower depths 
(effective gravity due to continuum radiation pressure). 

The work presented in this paper forms a basis for further 
development of dynamic photoionization models that will be 
covered in two following papers. In the first, the inclusion of 
chemistry and dissociation/formation processes allows a uni- 
fied treatment of the entire flow from cold, molecular gas, 
thr ough the PDR, and into the Hn region. Previous work 
by [Stoerzer & Hollenbach ( 1998 ) indicates that advection can 
have a significant effect on the position of the molecular hy- 
drogen dissociation front. In the second, divergent, transonic 
flows from strong-D fronts are modeled, which can explain 
the kinematic observations discussed in paragraph || above. 
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APPENDIX 



A. SIMPLIFIED ANALYTIC MODEL FOR WEAK-D IONIZATION FRONTS 

In this appendix, we develop a simple analytic model for a plane-parallel ionization front in order to explore the most important 
effects of advection upon the front structure. The principal simplification involved is the assumption that the gas temperature, T, 
follows a unique prescribed function of the hydrogen ionization fraction, x. It is a generic property of ionization fronts that the 
temperature tends to increase as one passes from the neutral to the ionized side of the front. In the analytic model, we assume 
that the increase is monotonic and has the form: 



T(x) 



1 +xfr 



(T m - T ), 



(Al) 



where T m is the limiting temperature in the ionized gas (x — > 1), Tq is the limiting temperature on the far neutral side (x — > 0), and 
/?t is a parameter controlling the sharpness of the transition. In reality, for moderate to high ionization parameters, the hardening 
of the radiation field as one approaches the ionization front causes the temperature to have a maximum for x somewhat less 
than unity (see, for example, Figure ^p). However, for weak-D fronts, equation ( |Al| ) is sufficient to capture the main effects of 
advection. 4 Also, it turns out that advection itself will modify the T(x) curve (see Section [|) but, again, we ignore this complication 
in the analytic model. We further simplify the model by only considering the ionization of hydrogen and ignoring the effects of 
radiation pressure, dust, and magnetic fields. 
With these approximations, the gas pressure is given by 



P = n{\ +x)kT, 



(A2) 



where n is the hydrogen number density. For a static front, the gas pressure will be constant, so it is a simple matter to calculate 
the variation of gas density, n, and electron density, n e , with ionization fraction: 



n(x) 



l+x \ r, 



T(x) 



n e (x) — xn{x). 



(A3) 



These are all plotted as solid lines in Figure 
j6t = 1/3 and To/T m = 0.02, which agrees within 



A18|. In this and all following examples, we use temperature-law parameters of 
0-20% with the temperature profiles of all the weak-D models in section f|. It 



4 For D-critical/strong-D fronts, on the other hand, the position of the temperature maximum is critical. 
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can be seen that, although the gas density declines with increasing x, the electron density is a monotonically increasing function 
of x (this is always true, whatever the value of /3j). 

In order to calculate the structure of a non-static, dynamic front it is necessary to consider the conservation of mass and 
momentum, which are given, in plane-parallel geometry, by 

nm H v = (Do and P + nm H v 2 — TIq, (A4) 

where v is the gas velocity and <t>o and n are the (constant) mass and momentum fluxes. It is convenient to define the dimen- 
sionless variable </>: 

<f> = £f~, (A5) 



which, by equations (|A4|), is related to the isothermal Mach number (vVt = v/c) via 



cf, = I (M + -L I , (A6) 



with the inverse relation 



M = 4> ± (<p 2 - l)' /2 , (A7) 



in which the - sign applies to a subsonic flow and the + sign to supersonic flow. 



Equation ( |A6| ) shows that cp > 1 everywhere and that 0=1 corresponds to the sonic point: M — 1 . A given ionization front 
can be characterized by the parameter <p m , which is the minimum value of <p anywhere in the front. This is achieved when the 
isothermal sound speed c attains its maximum value c m , and, in the case of the weak-D fronts considered here, corresponds to the 



maximum value of the (subsonic) Mach number: M m . For a temperature profile such as equation ( Al ), c m occurs at x = 1 and 
the sound speed varies with x as 

(T(x) 1 + jc\ 1/2 

c(x) = c m \ ^ 2 1 . (A8) 

It can be seen that all weak-D fronts with a given T(x) form a one-parameter family, characterized by their maximum Mach 
number M m . For a given M r . the corresponding <p m can be calculated from (|A6[) and then from ( A5) we have cf>(x) = cf> m c m /c(x), 



„ -in- — -— -- ^ — - -j.ij.j_r | j- in — - — — v| y — \ j \ / j " 1 ■ t ' i v 

which can be inserted into (A7) to give M(x), from whence we also have v(x) - M(x)c(x), «(x) = n m M m c m /v(x), and n e (x) 



xn(x). Henc e, the full structure of the ionization front as a function of x can be found algebraically. The results are plotted in 
Figure |Al8| for various weak-D fronts between M m — and 0.99. 



For M m - 0.2, the velocities are everywhere very subsonic and the density structure is hardly any different from the static case 
(lower panels). In this regime, the velocity rises approximately as u ^ c 2 Ai m /c m °c (1 + y)T(y). There is hence an initial brisk 
acceleration for y < 0.05, driven largely by the increase in T, followed by a slower, almost constant, acceleration, driven largely 
by the increase in ion fraction. For higher Ai m , the gas density contrast between the ionized and neutral sides increases and for 
M m > 0.7 this causes there to be a maximum in n e at an intermediate value of y. Solutions with M m > 0.7 also show a second 
episode of steep acceleration as y — > 1 . 

In order to apply these results, it is necessary to find the mapping between the ionization fraction, x, and physical position 
within the ionization front. For this, it is necessary to introduce more parameters than have so far been considered. The ionized 
gas is assumed to have an inner boundary z = 0, at which the ionizing flux is Fq, and with z being the distance into the ionized 
gas, measured from the illuminated face, in the direction of decreasing x. Globally, the flux of ionizing photons at z = must 
be balanced by the recombinations per unit area, integrated throughout the structure, plus the (constant) flux of hydrogen nuclei 
through the front: 

<h [*°° 

F = — + a(x)x 2 n 2 dz (A9) 



where a is the recombination coefficient, which we approximate as a power-law in the gas temperature: a = a m (T /T m )~ l . The 
global importance of advection on the ionization balance can be characterized by the relative magnitude of the two terms on the 
RHS of this equation, so we define a dimensionless "global advection parameter": 



Ad = (A10) 

F m H - O 



We also define a characteristic "Stromgren distance" as 



in terms of which, ( |A9| ) becomes 



zo = 1 + A f \ ax 2 n 2 dz, (All) 
a m n m Jo 



F = a m n 2 m zo. (A12) 



5 For R-type fronts, on the other hand, Mm is the minimum Mach number. For D-critical and strong-D fronts, A1 m = 1 and is no longer a turning point in At, 
not even a point of inflexion, although <fi m is still a maximum in (p. 
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At each point, z, within the structure we also have the local ionization equation: 



„2.„2 



d Oo dx 

F(l — x)no~ — — (xnv) — , 

dz m tl dz 



(A13) 



where F is the local value of the ionizing flux and <x is the mean photo-absorption cross-section for ionizing photons, evaluated 
by integrating over the ionizing spectrum at that point. The attenuation of the ionizing radiation is expressed by 



dF 

= —crn(l — x)F. 

dz 

Equations ( A13 ) and ( A14 ) can be reexpressed in dimensionless form as 

dx 1 
df 



<x(l -jc) 



and 



with dimensionless variables: 



-J = [T*cr«(l -xj] 1 . 
ar 



(A14) 



(A15) 



(A16) 



. _ z 

Zo' 



f,-ln^. 
Fo 



Equations (A15) and (416) also make use of two dimensionless parameters: 



t, = n m cr Zo, 



and 



C m O"0 



in terms of which the global advection parameter, /l a( j, can be expressed as 



^adWtn 



' fadWln 



(A17) 



(A18) 



(A19) 



and whose significance is explored more fully in section || These two differential equations can be integrated nume rically to find 
the full solution for the ionization front structure in physical space. Note that for the static case (M m = 0), ( A15 ) is undefined 
and we instead have simply that the term in square bracke ts is z ero. 



We now present sample results from solving equations (A.15) and ( A.16) using parameters appropriate to an Hn region illumi- 
nated by an O star. We approximate the reduction in photoabsorption cross-section due to the hardening of the ionizing radiation 

field as (T — [l + 0.2 (f + f 2 )] . This is a good fit to the exact result from assuming the ionizing spectrum of a 40 000 K black- 
body and cr(v) oc y~ 3 , from which we also obtain d~o = 0.505 cr(v ) 3 x 10~ 18 cm 2 . We also assume a m = 2.6 x 10~ 13 cm 4 s _1 , 
c m = 10 km s _1 , which give £ a d = 11.5. For the parameter t„, we adopt the values 30 and 3000, corresponding to a low and a 
high ionization parameter, respectively, with the second being more representative of typical Hn regions. The results are shown 
in Figure H. In each case, curves of the electron density and gas velocity are shown for models with (right to left) M m = 0.0, 0.2, 
0.5, 0.7, 079, 0.99. For r* = 30 the direct effects of advection are considerable, since = 0.64A1 m , so that for the higher Mach 
numbers a substantial fraction of the incident ionizing flux is consumed by ionization of new atoms, which pushes the ionization 
front to the left. 

It can also be seen that the advection decreases the width of the ionization front. In the static front, the width is determined by 
the mean free path of ionizing photons in the partially ionized gas, which gives 5z ~ (ner) , whereas in the advective fronts the 
ionization fraction at a given value of the ionizing flux is smaller, leading to a sharper front (4 times narrower in the near-critical 
case). Also, the electron density in the advective models has a peak at the ionization front, which is not seen in the static models. 
For t» = 30, we have A- d d = 0.0038vVt m , so the direct effects of advection on the global properties of the model should be very 
small. Nonetheless, the ionization front position varies by about 5% between the static model and the Ai m = 0.99 model, which 
is due to the effect of the electron density peak. 

In order to investigate the effects of advection on the emission-line properties of the nebula, we consider a generic recombina- 
tion line with emissivity 

Vrcc(.x) =A rec n e n j+1 T^, (A20) 
where n /+ ; is the number density of the recombining ion, and a generic collisionally excited line with emissivity 



A co \n e n j e 



-EjkT 



1 + B co in e T 



-1/2 



(A21) 



where E is the excitation energy of the upper level and Z? co i is the collisional de-excitation coefficient. The ion density can be 
assumed to be rij oc n e — xn for sin gly-io nized ions and rij oc (1 - x)n for neutral atoms. 



The results are shown in Figure A19 , which compares the line emissivities as a function of radius for the r« = 30 model at 
two values of 0.0 and 0.99. The total emission from all the lines is significantly reduced in the nearly D-critical model with 




z 

Fig. A19. — Line emissivity structure of static and dynamic ionization fronts (r» = 30), in which different line types corresponding to a generic recombination 
line (solid), optical collisional line of an ionized species (dashed), and of a neutral species (dotted). The set that peaks to the right is from a static model 
(M m = 0.0) while the set that peaks to the left is from a nearly D-critical model (M m = 0.99). 



respect to the static model due to the smaller depth of the ionized zone. On the other hand, the relative intensities integrated over 
the entire structure change by less than 10%. The peak of the neutral collisional line is much sharpened in the advective model 
due to the narrowing of the ionization front. In addition, the singly-ionized collisional line and the recombination line both show 
peaks in their emissivity at the ionization front, which are due to the electron density peak there (see Figure |J). 

The line emissivity can be combined with the velocity structure of the front to create synthetic emission line profiles. We 
assume that the front is observed face-on from the ionized side, so that the lines are all blue-shifted with respect to the neutral 
gas, which is assumed to be stationary. We calculate the emergent intensity profile, I(u), of an emission line ignoring any optical 
depths effects: 

I(u) = Tj(z) exp dz, 



2A 2 (z) 



(A22) 



where u is the observed velocity and A is the thermal Doppler width, calculated at each point in the structure assuming typical 
atomic weights of 16 for the collisional neu tral li ne, 14 for the collisional ionized line, and 1 for the recombination line. 

The resultant profiles are shown in Figure A20 for the collisional lines. 6 The increasing blue-shift of both lines with increasing 
maximum Mach number is evident. For the singly-ionized line, there is hardly any dependence of linewidth on the advection 
strength because most of the emission comes from the almost fully ionized gas, which shows only small velocity gradients. For 
the neutral line, on the other hand, the nearly critical model shows a much broader, double-humped line. The redder component 
of the line is due to the emissivity peak around x = 0.5 and is little changed between M m = 0.5 and M m = 0.99, whereas the 
bluer component in the M m = 0.99 model comes from the nearly fully ionized gas and is hence relatively stronger in the model 
with the higher ionization parameter (t» = 3000). 

The beha vior of the mean velocity, v, and RMS 7 velocity width, cr, of the collisional lines as a function of M m is shown in 
Figure A21. 



The recombination line has a very similar emissivity profile to the singly ionized collisional line and the low atomic weight only serves to smear out the 
details of the line profile. 

7 For a Gaussian line profile the full-width-half-maximum, Ao, is related to the RMS width by Av = 2.306<x. 
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Fig. A20. — Predicted line profiles from analytic ionization front model for {a) t, = 30 and (b) t, = 3000. Top panels shows a generic optical collisional line 
of a neutral species while bottom panels show the same for a singly-ionized species. Different line types correspond to M m = 0.0 (solid line), 0.5 (dashed line, 
and 0.99 (dotted line). 
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Fig. A21. — Predicted mean velocities and RMS widths (in km s ') of collisionally excited lines from the analytic ionization front model as a function of M, 
for (a) r, = 30 and (b) T. = 3000. 



B. TREATMENT OF IONIZATION LADDERS 
The rate of change of the fractional abundance of a particular ionization state is given by 



•i" j 



+2> 



(Bl) 



where the /?,_> ; - are the rates for ionization (where j is a higher state than ;) and recombination (where j is a lower state than i). G,- 
and S i cater for processes not included within the ionization ladder, and are respectively the source of ions from such processes 
and the sink rate into them. 
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In previous versions, Cloudy treated the ionization ladders in isolation, with G, = S ; = 0, and Z?y + only for processes 
coupling neighboring ionization states, 

[% ./-/- 1 

Ri-> j = \lj i = j + 1 (B2) 
1 otherwise 



In this case, equation ( [Bl| ) in equilibrium (d/dt 



-» 0) gives the equations 

= rc2??i - nil i 
= n 3 7? 2 -n 2 {I 2 +K{) 
= n 4 ^3 + n 2 I 2 - m(I 3 + %£) 
Q = n N -\I N -\ - nnHn-i 



(B3) 
(B4) 
(B5) 
(B6) 



for the abundances «, of an A^-state ionization ladder. These equations yield a simple expression for the relative abundances of 
neighboring ionization states, 

W«i = (B7) 

The overall abundance of each ionization level can be found using this relation together with a sum rule for the conserved total 
abundance of the species. 

This analysis does not apply if we require a time-dependent solution, or there are more complex interactions between levels 
(such as the Auger effect), or external sources and sinks of ions (resulting, for example, from molecular processes). 

If we assume a general form for all these additional terms, we are left with a computationally expensive NxN matrix problem to 
solve. However, the largest coefficients in the matrix derived from equation (Bl ) will either be the ionization and recombination 
rates, for which we know that a simple solution is possible, or possibly the time-dependent terms (as, for instance, in non- 
equilibrium cooling behind a shock). This suggests that we should be able to find a solution to the ladder equations efficiently 
using iterative techniques. 

We rewrite equation (Bl) as 

Y j A, i , l! -}jA i! -A i! ), li -b, (B8) 



We separate the matrix Ay into two parts, a tridiagonal component Ay and the remainder Ay. The components of Ay will in 
general be far smaller than those of Ay, so we can use the iterative scheme 



n" +1 =A~ X (b-A.n") 



+ A- 1 Co- 



A.n") 



(B9) 



to converge to the solution to the ionization state. In particular, in the nonlinear system we are treating, the coefficients in A and b 
will themselves be functions of the ionization state of the gas, and so it suffices to take a single step of the iterative scheme (B9) 
before these values are updated. 

There is one problem with this treatment. In the limit of small advection, A becomes singular. In this limit, the solution we 
require is the null eigenvector of A, and as in the previous treatment we can set its magnitude using an additional normalization 
constraint. However, the rounding error in the summation of ionization and recombination terms on the diagonal of A can lead 
to the numerical solution of equation (B9h having negative abundances for states which are substantially less abundant than their 
neighbors. The eas e with which the s olution (B7) is found suggests that this is not unavoidable. By re-writing the standard 
tridiagonal solver in ^ress et aT ( 1992 ) to treat matrices in the particular form of A (and providing the ionization, recombination 
and diagonal sink vectors to this revised algorithm without summation), a near-cancellation is avoided. The resulting scheme 
gives solutions which are manifestly positive, given the physical limits on the signs of the various vector elements. 

C. MAGNETIC FIELD 

Magnetic field effects can now be included in Cloudy simulations by specifying the magnetic field strength and geometry at 
the illuminated face. Both an ordered field and a 'tangled' field may be specified, although currently the ordered field is restricted 
to plane-parallel slab models with advection. The tangled field may be used in any geometry and with advective or static models. 

The tangled field is assumed to provide an isotropic magnetic pressure. In addition to the field strength at the illuminated 
face, Btangied.o, the effective magnetic adiabatic index, y mag , must also be specified. This determines the response of the field to 
compression of the gas: 



B 



tangled 



tangled,0 



(Cl) 



A value y mag = implies a constant magnetic field strength throughout the model, whereas y mag = 4/3 (the default) corresponds 
to conservation of magnetic flux and is what would be expected in the absence of dynamo action or magnetic reconnection. 

For the ordered field one must specify a component B z that is parallel to the integration direction through the slab and a 
component B, that is transverse to the integration direction. The parallel component B z is constant throughout the slab, while the 
transverse component is a function of the varying gas velocity, v: 
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where v* A is a characteristic speed ( Williams & Dyson 2001): 



4np v 

Magnetic pressure is included in the gas equation of state, having the form 



tangled D t "z i -2 

Anag = —3 + Z d y" e Cm ' 

57T 07T 



The magnetic contribution to the enthalpy density is given by 



